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SYMBOLS 


/ stream function 

S enthalpy 

u u(x,y)=f' F (x,y)-fX°°) 

v v(jc,jO = S F (x,y) - S(°°) 

w w{x,y)=f" p{x,y) 

X x that satisfies equations (7) through (11) 

x x=f"(o) 

Y y that satisfies equations (7) through (11) 

y y = S'(o ) 

z z(x,y) = S’ F {x,y) 

(3 pressure gradient parameter 

tj measure of distance from wall 

Subscripts 

F evaluation at a finite value of tj 

th 

i i = 1 , 2, or 3; i numerical value 

u partial differentiation with respect to u 

v partial differentiation with respect to v 

w partial differentiation with respect to w 

z partial differentiation with respect to z 

Superscripts 

( )' differentiation with respect to rj 

C) 


iterated value 
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SUMMARY 


A numerical method for the solution of large systems of nonlinear differential equations of the 
boundary-layer type is described. The method is a modification of the technique developed by 
Nachtsheim and Swigert for satisfying asymptotic boundary conditions. The present method 
employs inverse interpolation instead of the Newton method to adjust the initial conditions of the 
related initial-value problem. This eliminates the so-called perturbation equations required in 
Nachtsheim’s and Swigert’s application of Newton’s method. The elimination of the perturbation 
equations not only reduces the user’s preliminary work in the application of the method, but also 
reduces the number of time-consuming initial-value problems to be numerically solved at each 
iteration. Thus, although the rate of convergence of inverse interpolation is in general less than that 
obtained by the Newton method, the total number of required initial-value solutions may be 
reduced. For further ease of application, the solution of the overdetermined system for the 
unknown initial conditions is obtained automatically by applying Golub’s linear least-squares 
algorithm. The relative ease of application of the proposed numerical method increases directly as 
the order of the differential-equation system increases. Hence, the method is especially attractive 
for the solution of large-order systems. After the method is described, it is applied to a fifth-order 
problem from boundary-layer theory. 


INTRODUCTION 


Many problems in boundary-layer theory are posed mathematically as nonlinear, two-point 
boundary-value problems on an infinite interval. The so-called initial-value or “shooting” method, 
which does not require linearization of the equations, has been successfully applied to this class of 
problems. Basically, the shooting method adjusts the initial conditions so that the solution of the 
related initial-value problem satisfies the required boundary conditions at the outer boundary point. 

Of fundamental importance in the application of the shooting technique is the manner in 
which the adjustment of the initial conditions is made. Many adjustment schemes have been 
proposed (see ref. 1). In particular, Nachtsheim and Swigert (ref. 2) have developed an iterative 
adjustment scheme for boundary-layer problems with asymptotic boundary conditions. In that 
scheme, the differential equations are numerically integrated to some finite outer boundary point. 
The initial conditions are adjusted so that the errors between the computed values and the 
corresponding asymptotic values of both the functions with asymptotic boundary conditions and 



their derivatives are minimized in the least-squares sense. The iterative adjustments are continued 
until the solution converges within some prescribed accuracy. The precision of the solution is 
continually improved by progressively extending the range of integration. In this scheme, Newton’s 
method is applied to the nonlinear equations. However, Newton’s method requires the rate of 
change of the original differential equations with respect to each of the unknown initial conditions. 
Therefore, for each unknown initial condition, a concomitant differential-equation system (the 
so-called perturbation equations and their associated boundary conditions) must be formally 
constructed and numerically integrated simultaneously with the original differential-equation 
system. 

This paper describes a numerical method that modifies the Nachtsheim-Swigert adjustment 
scheme so that the perturbational differential-equation systems are eliminated. The numerical 
method is similar to the one proposed by Warner (ref. 3) for the finite interval. The elimination of 
the perturbation equations relieves the user of the burdensome task of formally generating and 
coding these equations for each problem, while simultaneously reducing the total number of 
numerical integrations to be performed. For additional convenience and accuracy, the linear 
least-squares solution for the unknown initial conditions is achieved using Golub’s method (ref. 4) 
of orthogonal transformations on the original system. So the task of finding the so-called normal 
equations, and then solving these linear algebraic equations for the least-squares solution, is also 
eliminated. The method is described and is applied to a specific problem in boundary-layer theory. 


NUMERICAL METHOD 


The method is best described by applying it to a specific problem. Consider Stewartson’s 
boundary-layer equations (for which Cohen and Reshotko, ref. 5, found similar solutions) with unit 
Prandtl number, and with pressure gradient parameter j3 = 1/2: 

/'"+//" + 1/2(5+ 1 -/' 2 ) = 0 (1) 

S"+fS' = 0 (2) 

f(o) = 0 , f(o) = 0 , S{o) = S w (3) 

/'(°°) = 1 , SM = o (4) 

where S w characterizes a particular solution. The generalization of the method to any order 
differential system follows naturally from the application to this fifth-order problem. 

For notational convenience, define the unknown initial conditions by 

x=f"(o ) (5) 

y = S\o ) (6) 

To apply the shooting method, initial guesses are made for x and y. These estimates, along with 
equations (1), (2), and (3), form an initial-value problem, which can be integrated by a numerical 
integration algorithm (e.g., Adams-Moulton) to some specified finite value of the independent 
variable p. Denote this outer boundary point of the finite integration range as r\p. Further, let a 
function subscripted with F denote the numerical value the function takes on when evaluated at r\p. 
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Define 


u(x,y) = f pix, y)~f\ °°) 
v(x,y) = Spix.y) - S(°°) 
w{x,y) =f"f(x,y) 
z(x,y) = Sj7(x,y) . 

The boundary-value problem of equations (1) through (4) could be solved by determining the 
unknown initial conditions x and y such that 


u(x,y) = 0 

(7) 

o 

ii 

. — . 

(8) 

w(x,y) = 0 

(9) 

z{x,y) = 0 

(10) 


Note that equations (7) and (8) require that the dependent variables f and S approach their 
boundary values in a finite domain, whereas equations (9) and (10) give the Nachtsheim-Swigert 
criteria, namely, the slopes of those variables also simultaneously approach zero. The solutions are 
not unique when only /' = 1 and S = 0 are specified at a finite r]p\ see reference 2. 

Instead of using a direct Newton iteration scheme to determine the roots x and y of the 
nonlinear system, equations (7) through (10) apply the inverse Aitken interpolation technique as 
generalized in reference 3. The functions u = u{x,y), v = v(x,y), w = wfx.y) and z = z(x,y ) may be 
inverted to give x = x(u,v), y = y(u,v) and also x = x(w,z), y = y{w,z), provided the Jacobians of the 
transformations do not vanish at any points in question. Denote the values of x and y when u, v, w, 

and z are zero as X and Y, and the values of derivatives dx/du , . . . , dy/dz when x = X and y = Y 

as X u , . . . , Y z . Expand x{u,v), y(u,v ) and x(w,z), y(w,z) into their corresponding Taylor series for 
two variables and then truncate after linear terms to obtain: 

x = X + uX u + vX v (11) 

y = Y + uY u + vY v (12) 

and 

x = X + wX w + zX z (13) 

y=Y + wY w +zY z (14) 

Because each of equations (11) through (14) contains three unknowns (e.g., X, X u , and X v in 
equation (1 1), it will be sufficient to specify three sets of initial trial values, (X! j-q ), (x 2 jq) an d 
(x 3 ,y 3 ). Equations (1) and (2) with conditions (3), (5), and (6) can be integrated to obtain (U]^,) 
and (wj.z,), (u 2 ,v 2 ) and (w 2 ,z 2 ), and (n 3 ,v 3 )and (w 3 ,z 3 ). Substituting these values into equations 
(11) through (14) and rearranging leads to the matrix equation 
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As in reference 2, the least-squares solution is the accepted solution of this overdetermined 
system of 12 equations with 10 unknowns. The usual method is not used— that is, determining the 
system of so-called normal equations and then solving the resulting linear algebraic system to obtain 
the least -squares solution of the overdetermined system, as in equation (15). 

As a replacement, the method developed by Golub (reference 4) is used to solve the linear 
least-squares problem. Golub’s algorithm is based on applying orthogonal Householder transforma- 
tions to the original overdetermined system. (See ref. 4 for a description of the method, as well as 
an error analysis.) This approach for obtaining the linear least-squares solution is numerically more 
stable than the solution via the normal equations. 

Once the values X and Y are determined from the least-squares solution of equation (15), they 
can be used in equations (5) and (6) along with equation (3) to again numerically integrate 
equations ( 1 ) and (2) to obtain new values ( 17, ,F, ) and (w { ,F, ). These values then replace the least 
accurate equations in equation (15). The resulting least-squares system is then solved again to 
determine the improved initial conditions, X and Y. This iterative procedure is continued until the 
errors between the computed boundary values and the specified asymptotic boundary conditions 
are within a prescribed tolerance, that is, until u, v, w, and z are all less than a very small prescribed 
number. 


NUMERICAL RESULTS 


A computer program was written using the proposed technique. The program was coded in 
FORTRAN IV for the time-sharing system (TSS/360) on the IBM 360/67 computer at Ames 
Research Center. As a specific example, the computer solved the problem given by equations (1) 
through (4) with the parameter S w equal to -0.2. 

The estimates for the unknown initial conditions ( x,y ) required to start the solution were (0.8, 
0.1), (0.8, 0.1001), and (0.8001, 0.1). The solution converged in 11 iterations so that the absolute 
error between the computed boundary conditions and the asymptotic boundary conditions was less 
than 1 X 10~ 9 . The initial conditions for the converged solution were f" = 0.86228190 and 
5' = 0.1062283, which agree with the values (five decimal places) given in reference 5. The total 
central processing unit (CPU) time was 6.85 seconds. 
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For comparison, the same problem was solved using the Nachtsheim-Swigert method. The only 
modifications to the computer program were those involving adjustment of the initial conditions. 
With the original guess of (0.8, 0.1) for the missing initial conditions, the computed values were 
within IX 10 -9 of their corresponding asymptotic values in seven iterations. However, the CPU 
time was 9.03 seconds. The increase in computer time was caused by the integration of three 
fifth-order systems at each iteration (namely, the original system and the perturbational system for 
each of the two unknown initial conditions) whereas, in the proposed method, only the original 
fifth-order system was integrated at each iteration. 


CONCLUSIONS 


For boundary-layer problems modeled by a large system of nonlinear differential equations 
with asymptotic boundary conditions (e.g., multicomponent diffusion flame problems), application 
of the shooting method with an adjustment scheme requiring perturbation equations is laborious if 
not prohibitive. The use of an inverse interpolatory adjustment procedure, which eliminates the 
perturbation equations, appears as an attractive alternative. The technique is easily programmed 
since only the original differential-equation system is used. And, although ease of application is the 
major advantage of this method, there may also be a savings in computer time, as was illustrated in 
the example considered. The suitability of the method appears to increase as the order of the 
system increases. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif. 94035, June 7, 1972 
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